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We study a system of dipolar molecules confined in a two-dimensional trap and subject to an 
optical square lattice. The repulsive long-range dipolar interaction D/r 3 favors an equilateral trian¬ 
gular arrangement of the molecules, which competes against the square symmetry of the underlying 
optical lattice with lattice constant b and amplitude V. We find the minimal-energy states at 
the commensurate density n = 1/b 2 and establish the complete square-to-triangular transformation 
pathway of the lattice with decreasing V involving period-doubled, solitonic, and distorted-triangular 
configurations. 


Competing structures and effects of commensuration 
appear in numerous physical systems. Prominent exam¬ 
ples are atoms on surfaces, e.g., Krypton on graphite [1], 
vortices in modulated superconducting films [2], in pe¬ 
riodic pinning arrays [3], and in a BEC subject to an 
optical lattice [4], flux quanta in Josephson junction ar¬ 
rays [5], or colloidal monolayers on periodic substrates 
[ ]. A new realization of this physics is accomplished 
by assembling cold dipolar molecules [7] (e.g., KRb [ 
or RbCs [9]) in a two-dimensional (2D) optical trap and 
stabilizing them with the help of a perpendicular electric 
field [1 ]. Adding a square optical lattice provides an 
effective substrate potential which competes against the 
triangular lattice arrangement favored by the long-range 
repulsive dipolar interaction. As a result, the system 
is expected to exhibit a variety of different configura¬ 
tions as a function of particle density and strength of the 
substrate potential. In this paper, we find the minimal- 
energy states at commensurate density in the absence of 
quantum and thermal fluctuations and thereby establish 
the complete transformation pathway from the square to 
the triangular lattice. Contrary to previous studies, the 
cold molecule system, besides being clean, can be contin¬ 
uously tuned through various configurations by changing 
system parameters such as particle density and substrate 
potential amplitude. Even more, modifying the orienta¬ 
tion or number of lasers, the symmetry of the optical 
lattice can be changed. 

In the simplest case, the transformation pathway be¬ 
tween lattices with different symmetries may involve a 
sequence of other uniform lattices. An interesting sit¬ 
uation arises when new topological objects show up in 
intermediate non-uniform phases. The original ‘misfit 
problem’ between a particle lattice with lattice constant 
a and a periodic substrate with incommensurate period¬ 
icity b ^ a has first been formulated in one dimension 
(ID); these studies [11, 12] have shown that the locked 
system at large potential V (with particle separation b) 
smoothly transforms into the free lattice (with separation 
a between particles) at V = 0 via a non-uniform soliton 
phase, with soliton cores approximating the free phase 
separating regions of locked phase. The commensurate- 


incommensurate transition in the 2D analogue has been 
addressed by Pokrovsky and Talapov [13-15]; within 
their ‘resonance approximation’, the problem reduces to 
a ID one and the system develops a secondary struc¬ 
ture in the form of a soliton-line array. Going beyond 
the resonance approximation, we find that the square-to- 
triangular transformation in the dipolar system involves 
three separate transitions related to the formation of a 
period-doubled zig-zag lattice as well as two instabilities 
towards non-uniform soliton phases, see Fig. 1. 



FIG. 1: Gibbs free energy of optimal states (thick line), tri¬ 
angular at V = 0, distorted and rotated triangular (gdt) at 
small V, solitonic and period-doubled (<? P d) at intermediate 
V, and square for V > Vb. Below the critical potential V/ 0 ’ 1 ’, 
the period-doubled phase smoothly transforms into the tri¬ 
angular lattice via two soliton transitions involving different 
soliton arrays. The dashed line extrapolates the energy <? p d of 
the period-doubled phase. Dotted lines are energies of rigid 
triangular (A), isosceles ([>), and square (□) configurations. 


We consider a 2D-confined molecular gas with dipolar 
interaction D/r 3 between the molecules, 
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subject to an optical (substrate) lattice 
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The particles with density n = 1/b 2 equal to commen¬ 
surate filling (one particle per minimum) are located 
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at positions r* with distances r,; ? = |r,; — r^ |; the sub¬ 
strate potential involves two modes qi = (g, 0 ) and 
q 2 = (0,g) with q = 2ir/b. Starting with the sys¬ 
tem’s energy for N particles trapped within the area A, 
E(A , N) = E' nt +E suh , our task is to minimize the Gibbs 
free energy per particle 

g(p) = G(A, N)/N = [E(A, N) + P A]/N, (3) 

where the thermodynamic limit with fixed density n = 
N/A is implied. We choose to work at fixed pressure 
p rather than fixed chemical potential, as this seems a 
better approximation to the experimental setup where 
molecules are confined to a trap. At V = 0 , the 
molecules arrange in an equilateral triangular lattice 
with a lattice constant a = (^/^Y^b > b and height 
h = (3/4) 1 / 2 a < b ; the resulting misfit parameter then 
is s = b/h — 1 ss 0.0746. Given the purely repulsive in¬ 
teraction between molecules, the density n is related to 
the pressure p = (3/2)ne A , with e A = e A 4 « 4.446 e D 
the interaction energy per particle in the triangular lat¬ 
tice and e D = D/b 3 the dipolar energy (the prefac¬ 
tor is conveniently calculated with an Ewald summation 
technique [16]). Upon switching on a small but finite 
potential V > 0, the rigid lattice assumes an energy 
g A (V) = e A (U) +p/n « 11.115 e D + U, increasing with 
amplitude U as each substrate mode contributes with an 
average V/2 to the energy. For a large potential U, the 
molecules arrange in a square lattice with lattice con¬ 
stant b < a and an energy g n « 11.186 e D independent 
of U as all particles occupy potential minima. Besides 
the triangular and fully locked square lattices, a third 
low-energy configuration [1 ] is that of an isosceles trian¬ 
gular lattice (below called the 66 -lattice) with base b and 
a height b locked to one substrate mode (we have to break 
the symmetry and choose the mode along x) with an en¬ 
ergy 9> {V) = 11.136 60 + 17/2. The above expressions for 
<? A (V), 9> (V), and g D already provide a reasonable ap¬ 
proximation to the energy g versus potential U diagram 
as illustrated in Fig. 1 (dotted lines). 

Next, we account for deviations U; of the particle co¬ 
ordinates r, : = R[ a4t + Uj from regular lattice positions 
R[ a4t . Expanding Eq. (1) in the displacement field Uj, 
the energy g = g\ att + 6g picks up a term 

( 4 ) 

i>3 

with the elastic matrix ^"(R 1 / 44 ) depending on the cho¬ 
sen lattice; the substrate potential contributes a second 
term 6 e sub to Sg, Sg = Sg lnt + 6 e sub . 

For a large substrate potential U, the substrate en¬ 
forces a square lattice with particle positions R la4t = R.P. 
Since the true configuration at F = 0 is the triangu¬ 
lar one, the square lattice becomes unstable when de¬ 
creasing U. The symmetry-breaking instability is to¬ 
wards a period-doubled zig-zag phase [l! ] and is con¬ 
veniently analyzed in Fourier space; the elastic matrix 


$°(k) exhibits negative eigenvalues, with the most nega¬ 
tive one (j>x = —3.958 e D n located at the A' point (n/b, 0) 
of the Brillouin zone and describing a shear distortion 
(alternatively, the symmetry breaking involves the point 
(0, 7 r/ 6 )). The contribution <5e sub shifts all eigenvalues by 
Vq 2 / 2, thus stabilizing the square lattice. The instability 
occurs when the first eigenvalue crosses zero at 

U p = -(2/g 2 )<^« 0.201 e D . (5) 


Decreasing U below Up, the system transforms to a 
period-doubled phase with two molecules per rectangular 
unit cell, see Fig. 1, with lattice vectors R™ c = (26,0) and 
R™ c = (0,6) and molecular positions c\ = (0,zti) and 
c 2 = ( 6 , u 2 ) therein. Inserting these coordinates into the 
functional Eq. (3) , we use the Poisson summation formula 
replacing the real-space sum along y by the reciprocal- 
space sum over £q to obtain the energy 
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£Ki[2nt(2i — 1)] 

8ire D ^2 -~ cos(q£6) ( 6 ) 

i,e> o 

+(U/2)[1 — cos(gcr) cos(g6/2)] + const. 


with a = (u\ + w 2 )/2 and S = (iq — u 2 ). The modified 
Bessel function K\(z) oc e~ z decays rapidly and we can 
limit the sum in Eq. (6) to the term i = i = 1. Mini¬ 
mizing g p d with respect to <5, we find that cos(g 6 / 2 ) = 
(U/8A) cos(qcr) with 2 A = g p — < 7 >(U = 0) « 0.0496 e D , 
and the energy reads 

U 2 

5 P d(cr) = 9> (V) - ^ cos 2 (gcr). (7) 

For the homogeneous period-doubled phase, <r = 0, i.e., 
the molecules displace symmetrically around the sub¬ 
strate minima along y , and 8 = (b/n) arccos(U/ 8 A). The 
condition 6 = 0 provides us with the critical potential 
Up = 8 A ss 0.198 e D \ this is close to the previous result 
(5), confirming that terms with i > 1 or l > 1 in Eq. ( 6 ) 
are indeed small. The order parameter approaches zero 
as 8 « ±(-\/2 6 / 7 r) a /1 — U/Up, while 8 = ±6/2 at U = 0 
describes the 66 -lattice with energy g > . The ± signs refer 
to the two possibilities to break the symmetry when dou¬ 
bling the period, leading to twin configurations with zig¬ 
zag structures shifted by 6 along x. The period-doubled 
phase then exists in four versions, with the zig-zag struc¬ 
ture manifest along x or y, each with a twin shifted by 6 . 
The energy y p d(U) of this phase resides below the energy 
g>{V) of the singly-locked isosceles phase, see Fig. 1. 

Next, we focus our interest to weak substrate poten¬ 
tials U. The particle coordinates then deviate from reg¬ 
ular triangular lattice positions, i.e., R latt = R/ and 
r i = Rf ± u, in Eq. (4). For very small U, one can 
expand the substrate potential to linear order in the dis¬ 
placement [19] and minimize the correction 8g in Fourier 
space. The force field involves the two modes q Q of the 
substrate potential, folded back to the first Brillouin cell 
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of the particle lattice, q a — n a Ki — m a K 2 = p a , with 
Ki, K 2 the reciprocal lattice vectors of the (triangular) 
particle lattice, n a , m a are appropriate integers, and we 
have included a minus sign in the definition of p Q for 
convenience. The minimal-energy configuration is found 
by rotating the triangular particle lattice with respect 
to the square substrate potential and relaxing the con¬ 
figuration in the force field. For a small misfit parame¬ 
ter s, one of the vectors p a passes near zero, generating 
a large deformation (and accordingly large energy gain) 
as the inverse elastic matrix [<f> D ] -1 (k —► 0) oc 1/k 2 is 
large at small k. Within the resonance approximation 
[14, 15], only the dominant term in the relaxation de¬ 
riving from the small misfit vector, say pi = Ki — qi, 
is included, while the small correction due to the other 
mode is dropped. Within this approximation, the op¬ 
timal value of the angle <p between the symmetry axes 
of the particle lattice and the substrate (see Fig. 1) is 
given by the same formula as derived by McTague and 
Novaco [19] for the accommodation of a triangular lat¬ 
tice on a substrate with the same (triangular) symmetry 
but with a different lattice constant, <p s ~ \/vs. Here, 
v = (k — //)/(« + n) is the Poisson ratio, with p, and k 
the shear and compression moduli. For the dipolar in¬ 
teraction oc R~ 3 , one has v = 9/11 [15] (k = 10/x and 
p/n = (3/8) e A ) and accordingly p = 3.86°. The energy 
(to leading order in s) of the distorted triangular phase 
reads 

V 2 n 

Sdt =5*00- 04^2“ (1 + m/k) (8) 

and we find a sinusoidal distortion held evolving along 
the direction z || p! enclosing an angle 9 = arctanydz k, 
42.13° with the substrate lattice, i.e., near the diagonal. 

With increasing V, this periodic distortion becomes 
large, of order 6, and turns into a soliton array as first 
described by Pokrovsky and Talapov [14, 15] within the 
same resonance approximation. Adopting a continuum 
elastic description and retaining the full anharmonic form 
of the substrate potential, they showed that the solution 
u = u g + Up minimizing the Gibbs free energy combines 
a global deformation u g with a periodic modulation u p 
that accounts for the soliton array. The global defor¬ 
mation u g involves a rotation and a uniform shear de¬ 
formation, smoothly transforming the rotated triangular 
lattice at V = 0 + into the isosceles lattice locked to the 
substrate along the x-axis at large V. In our case, this 
isosceles triangular lattice (below called the 66'-lattice) is 
characterized by a height b (along x), while, in the ab¬ 
sence of the second substrate mode, the base b ' ~ 1.0173 b 
along the y-axis can be found from minimizing the Gibbs 
free energy density g(p) at fixed height b and pressure p. 

The analysis of the soliton structure in Refs. 14, 15 
starts from the triangular lattice at small V and makes 
use of the associated isotropic elastic theory. Here, we fo¬ 
cus on the first soliton entry into the bb '-lattice upon de¬ 


creasing V; it then is more natural to calculate the energy 
of the deformation v (defined relative to the 66 '-lattice) 
using the elastic theory of the 66 '-lattice, g bb , (v) = g p + 
9 k +g^, with the linear term g p = ( 7 x+p)(d x v x ) + (y y + 
p ) ( d y V y ) driving the system towards the triangular phase 
and g K = K x (d x v x ) 2 /2 + K y (d y v y ) 2 /2 + K xy (d x v x )(d y v y ) 
and g y = p x {d y v x ) /2 T p y {d x v y ) /2 T p X y(dyVx)(dxVy) 
the usual (quadratic) elastic terms [21]; the coefficients 
are again calculated using Ewald techniques. In this for¬ 
mulation, the substrate energy assumes the simple form 
e sub = (W/2)[2 — cos (qv x )] with n! = 1 /bb' the particle 
density in the bb'- lattice. Aligning the rotated coordinate 
system (z,z±) with the misfit vector pi, the soliton dis¬ 
placement v(z) derives from a ID sine-Gordon equation. 

Using the isosceles elasticity, we find the Pokrovsky- 
Talapov (PT) soliton first entering the 66'-lattice at 
U C PT ~ 0.0417 e D \ the displacement held evolves along 
9 ss 45.05° (9 « 42.13° in the original analysis in Refs. 
14, 15 based on an isotropic elasticity, although see [20]) 
and shifts the particle lattice by d « (—6,0.70 6). With 
decreasing substrate amplitude V, the soliton density n BOl 
rapidly increases, n sol 6 oc 1/| ln(l — U/U C PT )|; the con¬ 
figuration with strongly overlapping solitons at small V 
then is equivalent to the rotated and distorted triangular 
phase obtained from perturbation theory. 

The soliton array obtained within the resonance ap¬ 
proximation transforms the 66'-lattice to the triangular 
one, while our goal here is to study the transformation of 
the particle system from square to triangular. The soli- 
tonic instability then should appear on the background 
of the period-doubled phase, which requires us to include 
the second harmonic of the substrate potential into our 
analysis. We expect the first soliton entry in the period- 
doubled phase to occur at small V where we can treat the 
period-doubled phase as an isosceles 66-lattice distorted 
by the relative shift d = b/2 — S of the two sublattices. 
Inside the soliton, the amplitude of this short-scale dis¬ 
tortion S = (6/7r) arcsin(U cos(gu y )/8A) is slaved to the 
center of mass coordinate v(R) replacing the scalar vari¬ 
able a introduced above. We then have to minimize the 
energy 

69 = Jjf d2R { 9 bb( v ) + ( 9 ) 

+ ^ [1_C0S ( 2 ^ y )]}’ 

where gf b is the elastic Gibbs free energy [21] density of 
the 66-lattice. While the resonance approximation admits 
only one low-energy soliton, the full problem with both 
substrate modes present allows for several line-defects 
shifting the lattice by d j : k = (— jb , kb/2) with j, k in¬ 
tegers. Promising candidates reminding about the PT 
soliton are the ( j, k) = (1, k) defects, but a simple Ansatz 
with the shift d 0 i = (0, 6/2) should be tried as well, since 
the particles merely have to overcome the weak effective 
potential oc U 2 /64A C b/2 along the y-direction, see 
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Eq. (9). All these line defects fall into two classes, the 
domain walls with j + k assuming odd values and tak¬ 
ing the period-doubled phase from one twin to the other, 
S — > —<5, and the genuine solitons with j + k even and 
the same twin on both sides, 5^5. 

The determination of the critical substrate potential 
for the (0,1) domain walls is straightforward, 

n V Ky + Hy cot 2 9 

and provides the maximal value V r c <0,1) ~ 0.0753 e D > V PT 
at 9 = 90°, see Fig. 2. The analysis for the (1, k) defects 
is more involved and the results depend strongly on the 
type of elasticity theory chosen for the calculation, telling 
us that corrections due to anharmonicities are large. 

For this reason, a reliable conclusion on the relevant 
scenario requires an numerically precise computation of 
the defects’ Gibbs free energies. Starting from a vari¬ 
ational Ansatz, we relax the particle configuration nu¬ 
merically for periodic arrays with large separations be¬ 
tween the defects. Summing up terms along the direc¬ 
tion perpendicular to z reduces the problem to a ID one, 
but restricts the possible angles 9 to those appertain¬ 
ing to small Miller indices. The results for the (0,1) 
domain wall (extrapolated to the thermodynamic limit) 
are shown in Fig. 2; they agree well with the analytic 
ones, although the largest V/ (0,1) « 0.07415 e D is assumed 
at a different angle 9 = 45°. While the flat dependence 
on angle renders the optimal orientation of the domain 
wall poorly defined, the data shows that the optimal de¬ 
fect does not align with a symmetry axis of the isosceles 
lattice. This result is quite unexpected, as such a sym¬ 
metry alignement is predicted by the analytic calculation 
neglecting anharmonicities and has often been considered 
as natural in the literature [22]. Our numerical results 
[23] for the (1 ,k) defects show that these would appear 
at much smaller values of V; in particular, the (second) 
best result V c (1,3) (0 = 45°) ~ 0.0544 e D is found for the 
(1, 3) soliton, while the k = 2 domain wall and k = 1 
soliton are even worse with I/, (1,2) (45°) ~ 0.0501 e D and 
V/ 1 - 1 ’^ 0 ) « 0.0382 e D . 

The proliferation of (0,1) domain walls washes out the 
y-harmonic and dilutes the particles along the j/-axis, 
thereby establishing the 6f/-lattice; the transformation 
to the (rotated) triangular lattice at V = 0 + then in¬ 
volves an additional PT solitonic transition at lower V 
which smoothly eliminates the x-harmonic. The ana¬ 
lytic result for V PT again can be improved with a nu¬ 
merical calculation and we find a maximal critical po¬ 
tential V PT {9 ss 44.5°) ~ 0.046 e D , see Fig. 2; at this 
value of the substrate potential, the domain-wall phase 
has approached the 66'-lattice to within « 10 %, as mea¬ 
sured by the ratio of amplitudes A p of the periodic de¬ 
formation v p generated by the (0,1) domain wall array, 
A p (V c PT )/A p (V^ 1> ) = 0.019/0.25 « 0.08. 



FIG. 2: Numerical results for the critical substrate potential 
V c for first soliton entry versus angle 9. Shown is the data 
for the (0,1) domain wall and for the PT soliton evaluated at 
selected angles defined by small Miller indices (m, n); dotted 
lines are guides to the eye. The flat form V r /°’ 1) (0) renders the 
angle 9 for the first (0,1) domain wall entry poorly defined. 
Thin lines are the analytic results following from a continuum- 
elastic description for an isosceles lattice. 


Depending on the specific situation at hand, alterna¬ 
tive scenarios can be realized. All of these have to respect 
that a phase transition establishing an array of identical 
solitons with shift vector d [e.g., (1, k ) solitons or domain 
walls] necessarily has to be followed by a further transi¬ 
tion at lower V; since the global distortion field in the 
soliton array is slaved to d, the rotated triangular phase 
at V = 0 + cannot be reached without the appearance 
of other defects. The completion of the transformation 
may then involve the formation of a network of crossing 
solitons. Furthermore, if the most favorable solitons have 
close critical potentials and intersect with a negative en¬ 
ergy, the two smooth transitions can merge into a single 
first-order one. 

To conclude, we discuss the prospects for an exper¬ 
imental realization and detection of these competing 
structures in a cold molecule system. In order to serve 
as a classical simulator, quantum fluctuations have to re¬ 
main small. While in usual cold atom systems the latter 
are limited by the optical lattice, here it is the long-range 
interaction between the molecules that bound the zero- 
point motion. In estimating the importance of quan¬ 
tum fluctuations, we have to compare the interaction en¬ 
ergy e D with the recoil energy e r = H 2 /mb 2 . Evaluating 
the quantum parameter tq = e D /e r ~ 15 Zd[D 2 ]/b[nm] 
(with Z denoting the molecular mass and D the De¬ 
bye unit) for favorable but reasonable parameter set¬ 
tings (Z ~ 100, b ~ 500 nm, \/~D ~ 5 D), we obtain 
tq ~ 10 2 . This is substantially larger than the critical 
value tq = r s { ss 18 [10] marking the transition to the su- 
perfluid state where quantum fluctuations dominate [24] . 
Hence molecular systems can serve as classical simula¬ 
tors, although some renormalization effects due to quan¬ 
tum fluctuations may occur. Furthermore, sufficiently 
large amplitudes V must be reached for the optical lat¬ 
tice; in a recent experiment [25], dipolar molecules have 
been localized in deep wells V ~ 10 2 e r , which should be 
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sufficient to reach the critical value Vp. Finally, a promis¬ 
ing way to identify the various structural phases is via 
their different dynamical response under an applied force 
field f, with the square and period-doubled phases char¬ 
acterized by symmetric and asymmetric (reduced along 
y) pinning, respectively. For practical purposes, the ex¬ 
ponentially weak pinning of solitons can be neglected; 
the force held f induces a drive f • dj t along z and the 
resulting soliton motion generates a mass how along d j ^ 
which allows to identify the two solitonic phases. 
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